arXivil 504.03647V 1 [cond-mat.mtrl-sci] 14 Apr 2015 


Preventing rapid energy loss from electron-hole pairs to phonons 

in graphene quantum dots 

J.P. Trinastic/ lek-Heng Chu/ and Hai-Ping Cheng*^ 

^Department of Physics and Quantum Theory Project, 

University of Florida, Gainesville, Florida, 32611, USA 
(Dated; April 15, 2015) 


PACS numbers: 73.21.La, 68.65.Hb, 78.67.Hc, 78.67.Wj 


1 


In semiconductors, photoexcited electrons and holes (carriers) initially occupy 
high-energy states, but quickly lose energy to phonons and relax to the band 
edge within a picosecond 1| . Increasing the lifetime of carriers in light-absorbing 
materials is necessary to improve open-circuit voltage in photovoltaics pj], charge 
separation in organic solar cells jsl, and charge transfer in photodetection de¬ 
vices j^. Here we demonstrate long lifetimes over one hundred picoseconds 
for electron-hole pairs in graphene quantum dots (GQDs) due to large transi¬ 
tion energies and weak coupling to excitonic states below the fundamental band 
gap. This possibility for a large transition energy to bound excitons is due to 
graphene’s poor screening, illustrating a unique mechanism in this QD to oc¬ 
cupy higher-energy states for long timescales. GQD edges can be terminated 
with either armchair or zigzag carbon patterns, and this edge structure changes 
excited state lifetimes by orders of magnitude. These results indicate nanoscale 
control of carrier lifetimes in optoelectronics. 

Semiconductor quantum dots (QDs) exhibit discrete transition energies between excited 
states due to electronic confinement at the dot boundary. In optoelectronic applications, 
a photoexcited carrier will transition to lower-energy excited states due to energy loss to 
lattice phonons, known as phonon-induced relaxation (slow relaxation: long state lifetime; 
fast relaxation: short lifetime). Energy must be conserved during this process such that the 


total energy of coupled phonons matches the transition energies around 0.1 eV 


B- 


This 

so 


B. 


transition energy often mismatches available phonon energies, typically tens of meV 
that multiple phonons must facilitate the transition to lower energies, leading to inefficient 
carrier relaxation that should create long excited state carrier lifetimes in QDs (’’phonon 
bottleneck”). Electronic coupling between involved states must also be weak, which d^ends 
on the wave functions and the momentum of phonon modes inducing the transition 7|. 

Understanding how to engineer phonon bottlenecks can improve optoelectronic devices 
1 a. however the promise of .o„g exeited earrier lifetimes beyoad tens of pieoseeonds 
(ps) has not been realized experimentally in QDs due to defects and surface ligands [8| 


lOj. These structural factors introduce additional states that reduce transition energies and 


increase electron-phonon coupling, leading to shortened carrier lifetimes in the picosecond 
range 8|. Thus, finding QDs with a phonon bottleneck and determining how structural 


changes control its timescale are priorities for improving optoelectronic performance. 
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:e a phonon bottleneck in GQDs 
ll| . The Coulomb interaction is 


In this Letter, we demonstrate a unique method to creat 
by taking advantage of graphene’s low dielectric constant 
not effectively screened in graphene and a strong attraction between excited electrons and 
holes creates bound states, known as excitons, with energies signihcantly below the funda¬ 
mental band gap. Originally described in 1931 by Frenkel, the exciton is a fundamental 
many-body excitation with neutral charge, consisting of an electron-hole pair typically lo¬ 
calized within one molecule or unit cell. Manipulating excitons has provided insights into 
charge transfer states at organic interfaces 12j, optical gain and laser emission 13|, and 


the excited state properties of low-dimensional materials 


14| . The low screening in GQDs 


should create a large transition energy between higher-energy states and excitons, requir¬ 
ing multiple phonons to induce a transition. This could lead to long excited state carrier 
lifetimes if other excited states have weak electronic coupling to excitons. GQDs are syn¬ 
thesized with chemical methods that minimize defects and allow for nanoscale structural 


control that we show can tune this electronic coupling. 


We model phonon-induced relaxation of electron-hole pairs in GQDs using a state-of-the 
art reduced density matrix method that rigorously accounts for excitons (see Methods). We 
hrst consider a 132-carbon-atom GQD with carbon ligands (C132-L) studied experimen¬ 


tally 1 ^. Figure [T](a) shows the general excitation and phonon-induced relaxation process 
predicted by our model results in Figure [I](b). To match transient absorption (TA) measure¬ 


ments Iq, we initially excite an electron-hole pair with a 3 eV photoexcitation (I), creating 


a ’’hot” electron-hole pair weakly bound by their Coulomb attraction and delocalized across 
the edges of the GQD. Over time, the pair loses energy to phonons (II), rapidly transition¬ 
ing downward through single-phonon processes over the hrst 0.1-1 ps. Within 10 ps, the 
electron-hole pair transitions to a state near 2.25 eV (S 3 in Figure [11(a)- (c)) that corresponds 
to the absorption peak of the molecule near the fundamental band gap (Figure [ 11 (c)). 

The S 3 state, however, does not correspond to the lowest excited state in GQDs, since 
poor screening leads to two excitonic states (Si and S 2 in Figure [Tl(c)) 0.30 eV below the 
absorption peak that match the experimental rise in absorption near 1.85 eV (red curves in 
Figure [11(c); Exp A has a small peak at 2.05 eV to be explored in future work). The large 
transition energy to the excitons necessitates a few-phonon process that drastically slows 
carrier relaxation, giving S 3 a long 130 ps lifetime before transitioning to S 1 /S 2 , localized 
near the center of the dot. This 100 ps timescale is 1-2 orders of magnitude longer than 
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QDs 


relaxation in bulk graphene [17| and an order of magnitude longer than lifetimes in other 


18| . The bound electron-hole pair stays in the excitonic states for two nanoseconds 


before nonradiatively relaxing to the ground state through a many-phonon process. Both 
the S 3 (130 ps) and Si (2.2 ns) lifetimes match extremely well with 76-185 ps and 1.5 ns 16 | 
lifetimes observed at similar energies in TA measurements. The 100 ps lifetime provides the 
potential to extract excited carriers from both the higher-energy excited states and excitonic 
states (Figure [T]( a), green region) in photovoltaic or other applications, as ex 


under 100 ps are possible to acceptor materials such as titanium oxide 


19 


xaction times 
To our 
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knowledge, these are the hrst ab initio results demonstrating slow carrier relaxation to 
excitons in a poor-screening material like graphene. 

Finding ways to modify state lifetimes would provide a path to tailor relaxation dynamics 
for various optoelectronic devices. Therefore, we next show how nanoscale structural changes 
can tune excited carrier lifetimes. As shown in Figure [21 we explore three structural modi- 
hcations: a) removing the carbon-chain ligands on the experimental C132-L GQD (C132), 
b) changing the carbon edge geometry terminating the GQD between armchair (AG) and 
zigzag (ZZ) patterns, and c) increasing the GQD size from 1-3 nanometers. For example, 
ZZ216 refers to a GQD with 216 carbon atoms and zigzag termination edges. 

Looking at the effect of ligands. Figures [T](b) and[T](d) plot relaxation dynamics for G132- 
L and G132, respectively. Lifetimes at high-energy states are similar, however the S 3 lifetime 
(rgg) is shortened by half in G132 (69 ps) compared to G132-L (130 ps). More dramatically, 
the Si lifetime (r^J decreases by an order of magnitude from 2.2 ns to 332 ps in G132, 
indicating that ligands delay nonradiative recombination to the ground state and could 
increase fluorescence yields. Previous work has indicated ligands expedite relaxation by 


introducing defects or hybridized states 


22 j, however we hnd that ligands increase lifetimes 


through physical mechanisms we explore below. 

Modifying carbon edge type (Figure [21(b)) and size (Figure [21(c)) reveals a dramatic ability 
to tune the phonon bottleneck to excitonic states. All GQDs demonstrate an absorption 
peak at S 3 and two excitonic states (Si and S 2 ) 0.23-0.53 eV below the fundamental band 
gap, with smaller GQDs exhibiting larger transition energies (see Tables SI-SII). Gomparing 
GQDs similar in size, we plot excited carrier energy as a function of time for AG114 and 
ZZ150 in Figures [31(a)- (b). In both cases, the excited carrier relaxes to the absorption peak 
(S 3 ) under a picosecond. However, the S 3 lifetime for AG114 (t^^^^^) dramatically reduces 
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to 0.4 ps, breaking the phonon bottleneck and providing a fast relaxation channel to excitons. 
Across all AC sizes (Figure |3l^d)), S 3 lifetimes are never more than 2 ps and higher-energy 
state lifetimes are within 10-50 fs {E/Eg > 1-5, where Eq is the hrst excitation energy). In 
contrast, the S 3 lifetime is two orders of magnitude larger for ZZ150 = 19.9 ps) and 

further lengthened by decreasing ZZ GQD size = 39.1 ps, Figure [21(c)). Increasing 

ZZ GQD size leads to a sharp decrease in S 3 lifetime to 3.1 ps for ZZ216 (Figure |3](e)). 

An intriguing effect of size on lifetimes occurs for high-energy states in ZZ GQDs (Figure 
|3](e)). Wider energy spacing at E/Eq > 1.7 (Figure SI) leads to increased lifetimes as a 
function of size, reaching picosecond lifetimes in ZZ216 (Figure |3](e)). Multiexciton genera¬ 
tion (MEG) can occur in states with E / Eg > 2 if phonon-induced relaxation is slow enough 
to allow high-energy carriers to Goulombically interact with another electron-hole pair and 
decay into multiple low-energy carriers. This process leads to multiple carriers per absorbed 
photon and increased photovoltaic efficiency. MEG occurs on timescales of 100 fs to 100 ps 
in GdSe 23| and 250 fs in PbSe or PbS 2^ QDs, competitive with the high-energy lifetimes 


seen in larger ZZ GQDs, encouraging future study of MEG in these systems. 

The above results suggest that the phonon bottleneck to excitons has a strong sensitivity 
to nanoscale changes in structure. To explore the physical mechanisms behind these differ¬ 
ences, we consider the factors influencing transition rates in our model (see Methods and 
Supplementary Materials). Phonon-induced transition rates are expressed as 25l. l26|: 


27r 

kafj “^|Fq./3| E{u!qiI^^, 


( 1 ) 


where is the traus.tion energy between states n and A Vis is the e.eetron.e eonpl.ng 
[7(], which measures wave function coupling via phonon momenta calculated using quantum 
molecular dynamics (QMD), and E{ua/ 3 ) is the Franck-Gondon-weighted density of states 
H- E{uJa/ 3 ) depends on the Huang-Rhys factor (Sap), which measures the coupling strength 
of a given mode to a transition. The function E{u) peaks at energies corresponding to modes 
with nonzero Sap- Examining Vap and Sap as a function of phonon frequency will provide 
physical insight into the carrier lifetimes discussed above. 

We hrst consider the effect of ligands on the above factors. Figures |l|(a)-(b) plot Vap for 
the ground state and lowest 10 excited states (indexed 0-10) in G132 and G132-L, respec¬ 
tively. Vap is smaller for G132-L across most matrix elements, indicating that the ligands 
suppress electronic coupling. In particular, ligands weaken coupling to excitons = 2.3 
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meV, = 1-7 meV) and to the ground state = 2.5 nieV, = 0.8 nieV), 

thus reducing kap and rationalizing the longer S 3 and Si lifetimes in C132-L compared to 
C132 (Figures [U^b) and[T](d)). The reduced likely occurs because ligands stabilize the 
dot from motion perpendicular to the graphene plane (Figure S2). During QMD, edge C 
atoms in C132 oscillate by several Angstroms whereas ligands in C132-L limit this motion. 
Such substantial nuclear motion in C132 increases electronic coupling between states. 

Both C132 and C132-L couple to low-frequency breathing modes as well as 1300 cm“^ 
(0.16 eV) and 1600 cm“^ (0.20 eV) defect and optical G modes typical in graphene (Figure 
[5](a)) [27| . The weak coupling of the S 3 -S 2 transition to these modes (S '32 < 0 . 01 ), combined 
with the 0.30 eV transition energy to excitons, leads to a small F{uai 3 ) and inefficient 
few-phonon relaxation process that results in the 100 ps carrier S 3 lifetimes and phonon 
bottleneck. The ligands in C132-L introduce additional phonon modes at 1000 cm“^ and 
3000 cm“^ , however they only weakly couple to low-energy transitions and therefore do 
not significantly impact lifetimes (Figure [5](a)). These results, combined with the reduced 
electronic coupling, provide a rationale for the 100 ps lifetimes seen experimentally Q. 

Examining both Vap and Sajs also explains different relaxation patterns in AC and ZZ 
GQDs. We focus on V 32 and V 31 , highlighted by the dashed blue line in Figures |l](a)-(f), 
as these couplings determine relaxation rates from the absorption peak (S 3 ) to excitons (Si 
and S 2 ). Both AC and ZZ GQDs exhibit strong electronic coupling to excitons (IG^ > 10 
meV). Whereas Vap in AC GQDs does not significantly change with size = 12.8 

meV, _ ii Q meV), coupling to excitons slowly decreases with size in ZZ GQDs 

= 9.7-21.0 meV). These values are significantly larger compared to C132 and C132-L 
and explain the shorter S 3 lifetimes, indicating that geometry and asymmetry play a crucial 
role in minimizing coupling and relaxation rates. 

The difference in relaxation rates between AC and ZZ GQDs arises due to large changes 
in exciton-phonon coupling. Transitions in AC GQDs couple more strongly {Sap > 0.01) to a 
range of phonon modes, including optical G and defect modes as well as hydrogen modes over 
3000 cm“^ (Figure[5](b)). Strong electronic couplings over 10 meV also lead to the fast 10-100 
fs lifetimes for E/Eg > 1-5 across AC GQD sizes (indices 7-10 in Figure 111(c)). In contrast, 
the ZZ edge leads to extremely weak electron-phonon coupling, as shown in Figure [5](c). For 
the S3-S2 transition, optical G and defect modes exhibit Sap < 0.01 even for the smallest dot. 
This weak coupling to phonons reduces E{uap) by orders of magnitude compared to the AC 


6 



GQDs, leading to 20-40 ps lifetimes in ZZ GQDs. This weak exciton-phonon coupling makes 
ZZ GQDs a promising material to generate long lifetimes if functionalization or geometries 
can be found to reduce the strong electronic coupling to excitons. 

In contrast to AG GQDs, the S 3 lifetime decreases with increasing ZZ GQD size due to 
the fact that the electron-phonon coupling is already so weak in the smallest ZZ GQD. As 
dot size increases, two competing effects determine kap- The transition energy decreases 
with size, requiring less phonons to induce the transition and increasing rates. However, 
electron-phonon coupling also decreases with size, which reduces Sap and Fiojap)- Whereas 
in AG GQDs the decrease of 00^2 with size is balanced by reduced Sap, in ZZ GQDs the effect 
of a smaller transition energy dominates and leads to shorter carrier lifetimes. In addition, 
0732 decreases to 0.23 eV in ZZ216, matching the 0.20 eV (1600 cm“Q optical G mode and 
creating a unique resonance condition that increases F{uap) by an order of magnitude, 
leading to a single-phonon process and short 3.1 ps S 3 lifetime. Further increasing ZZ 
GQD size would likely increase lifetimes because the resonance condition would no longer 
be satished and Vap and F{uap) would continue to decrease. 

Our report of long excited carrier lifetimes in GQDs suggest that large transition energies 
and weak electronic coupling to excitonic states provides a new class of phonon bottleneck 
exploited in low-screening materials. The bottleneck timescale is sensitive to the nanoscale 
structure of the GQD due to changes in electronic and electron-phonon coupling, providing 
a playground to tune state lifetimes by orders of magnitude for many applications. Asym¬ 
metric GQDs and ligands (G132-L) minimize electronic coupling whereas ZZ edges minimize 
exciton-phonon coupling, both of which should be of interest for hot carrier extraction in 
photovoltaics. Gonversely, the rapid relaxation in AG GQDs could be benehcial in light- 
emitting applications. Future work should explore functionalization and different GQD 
structures as a method to minimize coupling to excitons and maximize carrier lifetimes. 

Methods 


Computational details. Electronic excitations in GQ Ds were obtained by solving 
for the twenty lowest eigenvalues of the Gasida equation 28(] within linear response, time- 
dependent density functional theory (LR-TDDFT) 29!]. All calculations have been com¬ 
pleted using the Gaussian09 package 30j, employing the B3LYP hybrid exchange-correlation 
functional within the adiabatic approximation and the 6-31G(d) Gaussian basis set. The 
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exact exchange contribution included in the hybrid functional reproduces the correct asymp¬ 
totic behavior of the exchange potential that is necessary to accurately produce excitonic 
excitations below the electronic band gap in LR-TDDFT. Ground and excited state conhg- 
urations were relaxed such that all atomic forces are less than 0.01 eV/A. 

Quantum molecular dynamics (QMD) simulations were performed to produce input struc¬ 
tures for calculating the nonadiabatic coupling matrix elements (lA^) and room-temperature 
absorption spectra (see below). Structures were initially relaxed at OK and brought to 300K 
slowly through velocity rescaling. After equilibrating at 300K using the canonical NVT en¬ 
semble, a 1 picosecond (ps) simulation was performed within the NVE ensemble with a 0.5 
femtosecond (fs) time step. Structures from the QMD simulation were then used as input to 
calculate LR-TDDFT excitations and nonadiabatic coupling using the B3LYP functional. 

Absorption spectra. The zero-Kelvin absorption spectrum, for each GQD is 

calculated by weighting a delta function at each LR-TDDFT excitation energy, = hwoo, 
by its corresponding oscillator strength, /qq,: 


foJihuj - hwoa), (2) 

Q 

where /oa = |\ho) and ITq.) are the ground and excited state TDDFT 

wave functions, rUe is the electron mass, and (ToIrlTo) is the transition dipole moment. The 
delta function is modeled using a Gaussian function with 80 meV broadening. 

Forbidden transitions with near-zero oscillator strengths at OK , known as dark states, 
can be activated by thermal motion at room temperature. To better compare to experimen¬ 
tal absorption spectra, we model these finite temperature effects for the C132-L GQD by 
calculating and /qq, over multiple time steps of a 300K QMD simulation. The thermally 
averaged absorption spectra, A{u), is then 

1 ^ 

= ( 3 ) 

n=l a 


where N is the total number of time steps in the MD run, and and represent 
excitations and oscillator strengths at the nth time step. The thermally averaged absorption 
spectrum is used for all results shown. 

Relaxation dynamics. We combine a procedure to calculate nonadiabatic coupling 
within LR-TDDFT 


31 


32| with the reduced density matrix method [2^ to examine excited 


state population changes over time after an initial photoexcitation. To the best of our 







knowledge, this is the first development of a computationally efficient method that can be 
used to investigate systems with hundreds of atoms while treating excited states rigorously 
using LR-TDDFT, which allows for a correct description of excitons. 

Assuming a Markovian phonon bath and applying the secular approximation, the Redfield 


equations 


26j describe population changes in the reduced density matrix constructed 


from LR-TDDFT eigenstates as described by the transition rates ka/s (Eqn[T]). The electronic 
nonadiabatic coupling (14^) can be expressed as [^: 




ih 


Mu 


(vl/4V,|vI/^)(p), 


(4) 


where |Ta) is the auxiliary wave function constructed using the Casida ansatz [3l|, Vfe and 
Mk are the nuclear gradient operator and mass for the fcth atom, and (p) is the expectation 
value of the momentum operator. We evaluate this electronic coupling by running quantum 
molecular dynamics (QMD) simulations at 300K and calculating Vap between excitations a 
and (3 at each 0.5 fs time step using finite difference methods. 14^ typically converges within 
150 fs (Figure S3), highlighting the computational efficiency of the present method. 

Assuming a quantum harmonic oscillator phonon bath, FiujnFt) describes the density of 
phonon states coupled to a given nonadiabatic transition 26l |: 


1 roo 

FM = 

ZTTfl J —oo 


(5) 


where G{t) = Efc 5^/? - l)(n(cnfc) + 1) + - l)n(cnfc)]. 5^^ is the Huang-Rhys 

factor, a dimensionless quantity representing the coupling strength of the fcth normal mode 
to the transition, and n{uk) is the Bose-Einstein distribution for a phonon with energy huk 
at 300K. F{u}ai 3 ) can be expressed in terms of the nuclear wave function overlap and is 


connected to decoherence and the quantum Zeno effect 




To examine the excited population as a function of time, we consider a system bathed in 
light with frequency cj at f < 0. After the light is switched off at t = 0, a nonequilibrium 
carrier distribution will relax due to interactions with phonons as described above. The 
population change in aap is then defined as Ap(E, t) = o'aa{t)6{Ea — E) — rjeg, where Ea 
is the energy of the ath excitation and is the equilibrium population. Figures such as 
[T](b) and[3](a)-(c) plot Ari{E,t) after the initial photoexcitation. 

Please see the Supplementary Materials for more detailed explanation of the above 


9 








methodology. 
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FIG. 1. Summary of phonon-induced relaxation process and optical absorption for 132-carbon-atom 
graphene quantum dot (GQD). a) Diagram illustrating physical processes relevant to our model (see 
text for details), b) Model predictions of the energy as a function of time and state lifetimes (r) of 
an excited carrier in C132-L after 3 eV photoexcitation. Colors from blue to red indicate population 
changes from 0 to 1, where 1 represents the entirepopulation in one state. Experimental lifetimes 


taken from transient absorption measurements 


16| are included in 


absorption spectra of C132-L GQD comparing experimental (Exp A 


parentheses, c) Normalized 


Ifil ] and B 


2l|]) to calculated 


results (A(a;)). Vertical lines indicate the calculated excitation energies of the system (heights do 
not indicate a physical quantity). The black line indicates the HOMO-LUMO electronic band gap. 
Yellow and blue isosurfaces on the GQD structures indicate electron and hole charge densities of 
each excitation, emphasizing the localized nature of Si and S 2 (near-indentical charge densities) 
near the center of the GQD compared to higher excited states, d) Population relaxation dynamics 
for GQD without ligands (C132) similar to (b). 
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c) Dot Size 



ZZ54 



2.7 nm 


FIG. 2. Structural modifications to GQDs investigated to determine their effect on excited carrier 
lifetimes, a) 132-carbon-atom GQD with (C132-L) and without (C132) carbon-chain ligands used 
in experiment to prevent dot aggregation, b) GQDs with carbon atom edge termination following 
an armchair (AG) or zigzag (ZZ) pattern, c) GQDs of increasing size with diameters ranging from 
1 to 3 nanometers. GQDs are labeled by their edge type (AC or ZZ) and number of carbon atoms. 
For example, ZZ216 refers to a GQD with 216 carbon atoms and zigzag termination edges. 
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FIG. 3. Summary of relaxation dynamics for GQDs with armchair (AC) and zigzag (ZZ) termi¬ 
nation edges. All phonon-induced relaxation occurs after a 1.6 E/Eg photoexcitation, where Eq 
is the energy of the lowest excitation. Colors from blue to red indicate population changes from 
0 to 1, where 1 represents the entire population in one state at a given energy and time. Figures 
a) through c) plot the energy of the excited carrier as a function of time for AC114, ZZ150, and 
ZZ54, respectively, and list lifetimes (r) of important states. In all cases, S 3 corresponds to the 
absorption peak close to the electronic band gap and Si is an excitonic state, d) and e) plot carrier 
lifetimes (r) as a function of E/Eq and size for AC and ZZ GQDs, respectively. 
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FIG. 4. Electronic nonadiabatic coupling (V^/j) for GQDs with different structural modifications. 
(a)-(f) plot Vap matrix elements for the ground state (state index 0 ) and first ten excited states 
(state indices 1-10) for the C132, C132-L, AG114, ZZ54, ZZ150, and ZZ216 GQDs, respectively. 
The matrix elements corresponding to coupling between states at the absorption peak, S 3 , and 
excitonic states, S 2 (V 32 ) and Si (V 31 ), are highlighted by the dashed blue line. 
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FIG. 5. Huang-Rhys factors (Sa/s) as a function of phonon frequency, a) C132-L and C132, b) 
AC, and c) ZZ GQDs of varying sizes. The left and right columns refer to Huang-Rhys factors for 
modes coupling to the Si-Sq and S3-S2 transitions, respectively. 
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RELAXATION DYNAMICS METHODS 


Reduced Density Matrix Formalism 

To calculate transition rates between electronic excitations, we employ the reduced den¬ 
sity matrix formalism that leads to the famous Redfield equations. Comprehensive descrip¬ 
tions of the theory can be found in [l| and j^. The formalism provides a computationally 
efficient method for calculating dissipative dynamics of a system embedded in a large, macro¬ 
scopic environment. Here, we will consider electrons to be the system of interest coupled to 
phonons from the surrounding lattice acting as the environment, or reservoir. 

If the system and the reservoir make up a closed system, then we can write the total 
Hamiltonian as a sum of the system component, Hs, reservoir component. Hr, and the 
system-reservoir coupling component, V: 

= Hs + Hr + V. (1) 

The system and bath can each be described by their respective eigenvalue equations: 

( 2 ) 

HR\<Pi) = E, !</.,), (3) 

where |\['q) and Ea correspond to the system eigenstates and eigenvalues, and | 0 j) and Ei 
correspond to the reservoir eigenstates and eigenvalues. The density operator corresponding 
to Htot can be expressed as 

p(^) = (4) 

n 

where |<I)„) is a pure state of the closed, combined system and Wn is the corresponding 
probability. However, we are only concerned with the system degrees of freedom, so we can 
define a reduced density operator (RDO), d(t), by taking the trace over the reservoir space: 

a{t) = Trnlpit)] = ^ {(pi{t)\p{t)\<pi{t)). (5) 

i 

The RDO now only depends on the system’s degrees of freedom and allows us to understand 
its dynamics due to coupling to the reservoir. Moving to the interaction representation and 
assuming that the system-reservoir interaction is small enough that the system does not 
change the bath over time, the equation of motion for the RDO expanded to second order is 

= -jTrR[V\t), a\Q)pB\ ~ ^ dTTrB[V^{t), [V^{t - r), a^{t - t)pb]], (6) 

2 



where ps = represents the equilibrium thermal distribution of the 

reservoir (unchanged by the system) and [3 = l/fc^T, where ks is Boltzmann’s constant and 
T is the temperature. The hrst term on the righthand side of Eqn. E] involves a thermal 
average of V, which will only lead to an overall shift of the Hamiltonian and does not 
introduce dissipative effects. The second term contains information about the interaction 
between the system and reservoir that leads to irreversible energy dissipation from the 
system. Therefore, we will neglect the hrst term in the following discussion. 

We make two major assumptions to arrive at the hnal expression for relaxation dynam¬ 
ics between system eigenstates. First, the second term on the righthand side of Eqn. [6] 
introduces memory effects, seen by the fact that cr^(t) depends on a^{t — r). If we assume 
that the bath dynamics are much faster than the system dynamics, then u^{t) remains un¬ 
changed over the memory time scale and the effect can be ignored (r = 0), known as the 
Markov approximation. Applying this approximation and expressing the RDO in terms of 
the orthonormal system eigenstates, |Tq,), we arrive at the Redheld equations: 


dt 




( 7 ) 


where hwap = is the transition energy between eigenstates. Rap^iiv is a tetradic 

matrix that describes the rates of change of the elements of 

The second approximation, known as the secular approximation, reduces the number of 
elements of Rai 3 ,^u we need to calculate. Examining Eqn. [TJ we can see that the righthand 
side depends on an exponential term that oscillates rapidly. Over a long time average, the 
righthand side will vanish unless Uap — = 0. Imposing this constraint, the secular ap¬ 

proximation reduces Eqn. [7] to two cases: (1) a = /5,/i = z/, corresponding to population 
transfer between eigenstates, and {2) a ^ f3, a = p, f3 = u, corresponding to coherence de¬ 
phasing that describes phase differences between eigenstates. We can now write the Redheld 
equation for each of the two cases: 

1) Population relaxation: 


dt 




Raa,iifi dafj, 






( 8 ) 

(9) 


3 




2 ) Coherence dephasing: 


= -Rc.p,c.p(Tap{t) 


Raf3,af3 ^ ^ r) (^«At ^/3 m) "h ^O) 


where 70 is a pure dephasing constant (Sj]. The rates ka /3 can be shown to be 




27r 


kap = — y^Ji\{(j)i'^a\V\(j)j'^0)\‘^6{hWafi - hWij] 




( 10 ) 

( 11 ) 


( 12 ) 


= X] l(‘)’i'i'i>ir|')’j'i'/>>l7ii^os - fiw.j)). (13) 

j 

where /j = and, in the second line, we have rewritten the initial distribu¬ 

tion of quantum reservoir states as a thermal average over initial conditions, = Si/*- 
The sum over hnal reservoir states preserves conservation of energy. The principle of detailed 
balance ensures probability conservation and provides the following relationship between for¬ 
ward and backward transition rates: 


kafi = (14) 

In the present study, we are only concerned with population transfer between eigenstates of 
the system Hamiltonian, and therefore we only solve Eqn[H]and ignore coherence dephasing. 

In the next section, we describe the system-reservoir interaction V that leads to nona- 
diabatic coupling between system eigenstates and determine how to calculate the matrix 
elements in Eqn. [121 
Nonadiabatic Coupling 

Within the Born-Oppenheimer approximation, the total wave function is expressed as 
a product of electronic and nuclear components and the nuclear kinetic energy operator 
is neglected. This leads to an electronic Hamiltonian that includes the electronic kinetic 
energy, electron-electron repulsive potential, and the electronic-nuclear attractive potential. 
The eigenstates of the electronic Hamiltonian depend only parametrically on the nuclear 
coordinates. Transitions between these adiabatic eigenstates can be described by including a 
nonadiabatic coupling perturbation dehned by the nuclear kinetic energy operator, modeling 
transitions induced by lattice vibrations. The matrix elements of this nonadiabatic coupling 
operator are 

7- ^ 



4 







where Mk is the nuclear mass and k runs over the K nuclear degrees of freedom, |0j) = 
\(f)])\(f)j)...\(f)f). If we only include terms to hrst order in the nuclear gradient operator, Eqn. 
becomes 

Vaf},ij = ^ —(d'„|Vfc|d'/3)(0i|VA:|0j), (16) 

k ^ 

such that the operator becomes a product of electronic and nuclear coupling elements. Pre¬ 
vious work has shown that Eqn. [16] can be rewritten to good approximation as 

ih 

n('^i Id). (17) 

k ^ k 


where (p) = —ih{(l)i\Vk\4’i) is the expectation value of the nuclear momentum operator, 
which follows classical equations of motion according to Ehrenfest’s theorem. Using this ex¬ 
pression to represent the matrix elements of the system-reservoir coupling V, our relaxation 
rate from Eqn. [13] becomes 

-T 7 -('i'»iwi>i's)(p)r E n (18) 

k ^ j k 


where the electronic nonadiabatic coupling is 

order of magnitude difference in electron and nuclear relaxation times, we apply a decorrela¬ 
tion assumption to separate the thermal averaging of the electronic and nuclear components 


B 


• Ig: 


kcs = - fa.,), (19) 

i,j k 

Thus, calculation of the transition rates between excitations reduces to the calculation 
of two major quantities: 1) the thermally averaged square of the electronic nonadiabatic 
coupling, and 2) the Franck-Condon-weighted phonon density of states (FC-DOS) 

, F{uap) = Yhi j fi rifc l(0f f^e Same term that appears in the absorption 

spectra but with a; = 0. Previous work has calculated transitions rates using Eqn. [19] but 

ignored the quantum nature of the nuclear wave functions, essentially setting the nuclear 

overlap to one 0.0: 

kafi = y^JiSjhWaiB - hjij). ( 20 ) 

id 


This approximation has given results matching experiment in some cases, however likely 
overestimates relaxation rates because it weights all phonon modes equally. Especially in 


5 









systems with a large number of bath modes and weak electron-phonon coupling, it is impor¬ 
tant to account for the nuclear overlap to only include modes coupled to a given transition. 
Nonadiabatic Coupling within LR-TDDFT 

We calculate the electronic nonadiabatic coupling matrix elements usingthe auxiliary 
many-body wave function originally proposed by Casida {q] and Tavernelli 1^ to describe 
excited state wave functions within LR-TDDFT. In this case, wave functions are expressed as 
linear combinations of singly excited Slater determinants of Kohn-Sham orbitals, determined 
from the solution of the Casida equation: 


*.) = yjxr + = T(A'r+ 


ma 


ma 


( 21 ) 


where m and a index occupied and unoccupied Kohn-Sham orbitals, respectively, and 
Qm are creation and annihilation operators, and are the normalized coefficients 
corresponding to an excitation from the mth to ath orbital and de-excitation from the ath 
to mth orbital, respectively, and |\ko) is the ground state Slater determinant. This auxiliary 
wave function has been shown to give the exact hrst-order nonadiabatic coupling between 
the ground state and any excited state and is an accurate approximation to the hrst- 
order coupling between excited states [lit]. The Vap matrix elements can be calculated by 
converting the nuclear gradient operator to a time derivative and using the hnite difference 
method between adjacent timesteps in a quantum molecular dynamics simulation: 

ih 




M, 




( 22 ) 


d 


ih 


2At 


+ At)) - (T„(t + At)|T;3(t))], 


where At is the 0.5 fs time step of the molecular dynamics simulation. At each time step, 
we calculate the nonadiabatic coupling between the many-body wave functions using Eqns. 
fl2T]) and fl22l) . In particular, we can decompose the many-body wave function overlap from 
the hnite-difference method in Eqn. (|2^ as a sum over singly excited Slater determinant 
overlaps: 

(T„(t)|T^(t + At)) = (23) 

yixrit) + Yrmxfit + ao+ y/(« + + a«)). 

ma nb 
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The singly excited Slater determinant overlap in Eqn. fl23|) can be decomposed using Slater- 
Condon rules and calculated using the overlaps between single-particle Kohn-Sham orbitals 
making up each determinant. These overlaps are nonzero when either the determinants are 
the same or differ by only one single-particle orbital. Thus, the nonadiabatic coupling matrix 
elements (V^/j) are found by weighting the Kohn-Sham orbital overlaps by the corresponding 
LR-TDDFT excitation coefficients and and summing over all single-particle 

transitions that make up each many-body excitation. 

The thermally averaged electronic nonadiabatic coupling that enters into Eqn. flT^ is 
calculated by averaging the value of Vk /3 calculated between adjacent time steps across a 
300K QMD trajectory consisting of N time steps: 

1 ^ 

= (24) 

i=l 

Pranck-Condon-Weighted Density of States 

We calculate the Franck-Condon-weighted density of states assuming that the phonon 
reservoir can be described as a harmonic bath and that the ground state normal modes are an 
accurate approximation for all excited state modes. In this case, nuclear wave functions can 
be expressed as eigenstates of a quantum harmonic oscillator and Fi^ojap) can be rewritten 

as |i|: 

Fiu^a) = (25) 

27rh J_^ 

and 


G(i) = - l){nM + 1) + (e”‘‘ - l)n(u)»)|, (26) 

k 

where cjk is the frequency of the fcth mode, t is time, and n{ujk) = — 1)“^ is the 

Bose-Einstein distribution that provides temperature dependence. is the Huang-Rhys 
factor, defined as 


Qk _ ,2 


(27) 


where fik is the reduced mass of the fcth normal mode and dk is the projection of the mass- 
weighted, configurational difference between excitations a and /3 onto the fcth mass-weighted, 
normal mode eigenvector. The Huang-Rhys factor measures the coupling strength between a 
given electronic transition and the kth mode and is the primary contribution to altering the 
FC-DOS compared to the complete phonon DOS. Modes that couple strongly with a given 
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TABLE I. Vertical excitation energies {Eq^) and oscillator strengths (/oa) for low-lying states at 
OK for armchair (AC) GQDs. 


AC42 AC114 AC222 


Excitation E 

Ta (eV) 

foa 

Eg- (eV) 

foa 

Eg- (eV) 

foa 

So ^ Si 

2.894 

0.000 

1.900 

0.000 

1.418 

0.000 

So^Ss 

3.028 

0.000 

1.970 

0.000 

1.462 

0.000 

So^Sg 

3.291 

0.000 

2.306 

1.328 

1.738 

2.029 

So —)• S4 

3.476 

0.722 

2.309 

0.000 

1.783 

0.000 

So^Ss 

3.624 

0.000 

2.462 

0.000 

1.873 

0.000 

So So 

3.834 

0.000 

2.739 

0.000 

2.139 

0.000 


transition corresponds to a larger peak at that mode energy in the FC-DOS. Multiphonon 
processes are also included in Eqn. [25] that result in diminished peaks at multiples of a 
given phonon energy. Finally, the reorganization energy, A, for a given electronic transition 
is defined as 

(28) 

k 

which measures the energy dissipated in the system to relax to its minimum after a vertical 
electronic transition. Small reorganization energies correspond to transitions between states 
with similar geometries. 

In summary, the calculation of phonon-assisted relaxation between TDDFT excitations 
requires three ingredients: 1) the reduced density matrix constructed from the eigenstates 
of the Casida equation within linear-response TDDFT; 2) thermally averaged nonadiabatic 
couplings Vap calculated along a QMD trajectory within the LR-TDDFT framework; and 
3) the Franck-Condon-weighted density of states (FC-DOS) calculated within the harmonic 
approximation and requiring the ground state normal modes and excited state configurations 
of each system. 




TABLE II. Vertical excitation energies {Eq^) and oscillator strengths (/oa) for low-lying states at 
OK for zigzag (ZZ) GQDs. 

ZZ54 ZZ96 ZZ150 ZZ216 


Excitation E 

(eV) 

foa 

Ka (eV) 

foa 

Ka (eV) 

foa 

Ka (eV) 

foa 

So^Si 

2.239 

0.000 

1.667 

0.000 

1.280 

0.000 

0.992 

0.000 

So^S2 

2.391 

0.000 

1.806 

0.000 

1.422 

0.000 

1.139 

0.000 

So^Ss 

2.905 

0.929 

2.179 

1.143 

1.688 

1.279 

1.324 

1.311 

So^S4 

3.122 

0.000 

2.367 

0.000 

1.818 

0.000 

1.379 

0.000 

So^Ss 

3.140 

0.000 

2.386 

0.000 

1.865 

0.000 

1.474 

0.000 

So^Se 

3.205 

0.000 

2.429 

0.000 

1.905 

0.000 

1.497 

0.000 


TABLE III. Transition energy {oJap), electronic nonadiabatic coupling (V^/s); and Franck-Condon- 
weighted density of states (F^ujais)) for sample transitions between LR-TDDFT excitations in 
C132-L and C132 GQDs. 


Transition (eV) 


Vai3 (meV) 


F(i 0 ^f,) (1/eV) 



C132-L 

C132 

C132-L 

C132 

C132-L 

C132 

Si-So 

1.918 

1.920 

0.794 

2.550 

0.073 

0.052 

S 2 -S 1 

0.060 

0.061 

20.075 

36.606 

0.139 

0.044 

S 3 -S 2 

0.295 

0.302 

1.649 

2.324 

0.134 

0.065 

S 3 -S 1 

0.355 

0.363 

1.729 

5.206 

0.098 

0.027 

S 4 -S 3 

0.069 

0.062 

16.733 

8.087 

0.270 

0.168 

S 4 -S 2 

0.364 

0.364 

1.691 

3.701 

0.194 

0.089 

S 4 -S 1 

0.424 

0.425 

1.905 

7.183 

0.110 

0.034 

S 9 -S 8 

0.071 

0.076 

2.156 

9.105 

2.198 

1.617 
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FIG. 1. LR-TDDFT excitations for various GQD types and sizes, a) Excitation levels for armchair 
(AG) GQDs. b) Excitation levels for zigzag (ZZ) GQDs. Excitation energy is given as a function 
of E/Eg, where Eq is the lowest LR-TDDFT excitation energy. Blue horizontal lines indicate 
excitation energies. The red excitation indicates a state of interest that changes its position relative 
to other states as a function of size in AC GQDs. C stands fo the continuum level, at which point 
transition energies between states are only several meV. 
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FIG. 2. Configuration of a) C132-L and b) C132 GQD at representative timestep in the 300K 
molecular dynamics trajectory. Medium gray spheres are carbon atoms in the body of the GQD, 
medium blue spheres are carbon atoms in the ligands, and small pink spheres are hydrogen atoms. 
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FIG. 3. a) Temperature and b) average square of the nonadiabatic coupling (< |V'a/3p >) as a func¬ 
tion of time from the 300K ab initio molecular dynamics simulation run within the microcanonical 
(NVE) ensemble for the C132-L and C132 GQDs. 
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